{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.4.1"
    },
    "colab": {
      "name": "406b_Problem_Sheet.ipynb",
      "provenance": [],
      "include_colab_link": true
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/john-s-butler-dit/Numerical-Analysis-Python/blob/master/Chapter%2004%20-%20Multistep%20Methods/4_Problem_Sheet/406b_Problem_Sheet.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "v9EAcof3thEm"
      },
      "source": [
        "## Problem Sheet Question 2a\n",
        "\n",
        "The general form of the population growth differential equation\n",
        "\\begin{equation} y^{'}=ty^3-y, \\ \\ (0 \\leq t \\leq 2) \\end{equation}\n",
        "with the initial condition\n",
        "\\begin{equation}y(0)=1.\\end{equation}\n",
        "For N=4\n",
        "\\begin{equation} y(x_1)= 0.5.\\end{equation}\n",
        "\n",
        "### 2-step Adams Bashforth\n",
        "\n",
        "The 2-step Adams Bashforth difference equation is\n",
        "\\begin{equation}w^{0}_{i+1} = w_{i} + \\frac{h}{2}(3f(t_i,w_i)-f(t_{i-1},w_{i-1})) \\end{equation}\n",
        "\\begin{equation}w^{0}_{i+1} = w_{i} + \\frac{h}{2}(3(t_iw_i^3-w_i)-(t_{i-1}w_{i-1}^3-w_{i-1})) \\end{equation}\n",
        "### 3-step Adams Moulton\n",
        "\\begin{equation}w^{1}_{i+1} = w_{i} + \\frac{h}{12}(5f(t_{i+1},w^{0}_{i+1})+8f(t_{i},w_{i})-f(t_{i-1},w_{i-1})) \\end{equation}\n",
        "\n",
        "\\begin{equation} w^{1}_{i+1} = w_{i} + \\frac{h}{12}(5(t_{i+1}(w^0_{i+1})^3-w^0_{i+1})+8(t_{i}w_{i}^3-w_{i})-(t_{i-1}w_{i-1}^3-w_{i-1})). \\end{equation}\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "collapsed": true,
        "id": "WuMO5Vh2thEp"
      },
      "source": [
        "import numpy as np\n",
        "import math \n",
        "\n",
        "%matplotlib inline\n",
        "import matplotlib.pyplot as plt # side-stepping mpl backend\n",
        "import matplotlib.gridspec as gridspec # subplots\n",
        "import warnings\n",
        "\n",
        "warnings.filterwarnings(\"ignore\")\n"
      ],
      "execution_count": 1,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "collapsed": true,
        "id": "fV7hT4BlthEr"
      },
      "source": [
        "def myfun_ty(t,y):\n",
        "    return y*y*y*t-y\n",
        "\n",
        "\n",
        "\n",
        "#PLOTS\n",
        "def Adams_Bashforth_Predictor_Corrector(N,IC):\n",
        "\n",
        "    x_end=2\n",
        "    x_start=0\n",
        "    INTITIAL_CONDITION=IC\n",
        "    h=x_end/(N)\n",
        "    N=N+2;\n",
        "    t=np.zeros(N)\n",
        "    w_predictor=np.zeros(N)\n",
        "    w_corrector=np.zeros(N)\n",
        "   \n",
        "    Analytic_Solution=np.zeros(N)\n",
        "    k=0\n",
        "    w_predictor[0]=INTITIAL_CONDITION\n",
        "    w_corrector[0]=INTITIAL_CONDITION\n",
        "    Analytic_Solution[0]=INTITIAL_CONDITION\n",
        "    t[0]=x_start\n",
        "    t[1]=x_start+1*h\n",
        "    t[2]=x_start+2*h\n",
        "    w_predictor[1]=0.5\n",
        "    w_corrector[1]=0.5\n",
        "   \n",
        "    for k in range (2,N-1):\n",
        "        w_predictor[k+1]=w_corrector[k]+h/2.0*(3*myfun_ty(t[k],w_corrector[k])-myfun_ty(t[k-1],w_corrector[k-1]))\n",
        "        w_corrector[k+1]=w_corrector[k]+h/12.0*(5*myfun_ty(t[k+1],w_predictor[k+1])+8*myfun_ty(t[k],w_corrector[k])-myfun_ty(t[k-1],w_corrector[k-1]))\n",
        "        t[k+1]=t[k]+h\n",
        "   \n",
        "    fig = plt.figure(figsize=(10,4))\n",
        "    # --- left hand plot\n",
        "    ax = fig.add_subplot(1,2,1)\n",
        "    plt.plot(t,w_predictor,color='red')\n",
        "    #ax.legend(loc='best')\n",
        "    plt.title('Predictor h=%s'%(h))\n",
        "\n",
        "    # --- right hand plot\n",
        "    ax = fig.add_subplot(1,2,2)\n",
        "    plt.plot(t,w_corrector,color='blue')\n",
        "    plt.title('Corrector')\n",
        "\n",
        "    # --- titled , explanatory text and save\n",
        "    fig.suptitle(r\"$y'=ty^3-y$\", fontsize=20)\n",
        "    plt.tight_layout()\n",
        "    plt.subplots_adjust(top=0.85)    \n",
        "    print('time')\n",
        "    print(t)\n",
        "    print('Predictor')\n",
        "    print(w_predictor)\n",
        "    print('Corrector')\n",
        "    print(w_corrector)"
      ],
      "execution_count": 2,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "scrolled": true,
        "id": "9eUbqJYjthEs",
        "outputId": "9bfd2efd-a628-4e0f-e68f-0c7b8e76bcdf",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 470
        }
      },
      "source": [
        "Adams_Bashforth_Predictor_Corrector(10,1)"
      ],
      "execution_count": 3,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "time\n",
            "[0.  0.2 0.4 0.6 0.8 1.  1.2 1.4 1.6 1.8 2.  2.2]\n",
            "Predictor\n",
            "[1.00000000e+00 5.00000000e-01 0.00000000e+00 4.75000000e-02\n",
            " 2.77084450e-03 2.63559724e-03 2.15353298e-03 1.76273674e-03\n",
            " 1.44282559e-03 1.18097388e-03 9.66644343e-04 7.91212442e-04]\n",
            "Corrector\n",
            "[1.00000000e+00 5.00000000e-01 0.00000000e+00 3.95833333e-03\n",
            " 3.19965681e-03 2.61937789e-03 2.14399600e-03 1.75489271e-03\n",
            " 1.43640563e-03 1.17571911e-03 9.62343270e-04 7.87691972e-04]\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAEdCAYAAAARsJF3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3debydZX3v/c8vM5hIkASFDASRUZm3YQgkkOz7HLUWPMcOaAdpPfV4LLbW1tZWai2Ptn30sfq0clSsVus5Dkg9Ptjiy5IECCQBEuYZw5hEhpDQMGbkev641m5WNntnT2uva629Pu/Xa73WdO+1futm57e/XOu+7itSSkiSJEnKxpUuQJIkSWolBmRJkiSpjgFZkiRJqmNAliRJkuoYkCVJkqQ6BmRJkiSpjgFZkiRJqmNAliRJkuoYkCWpjUTE70bEnRHxXO2yOiJ+oXRdkjSWhCvpSVL7iIjzgR3Az8iDHO8D/hg4NaV0Z8naJGmsMCBL6igRsQq4OqX0F6VraZSI2AL8aUrpq6VrkaSxwEMsJHWMiDgKOAPYUOC9/yAiUkS8t4GvOT4iLgCmAqsa9bqS1OkmlC5AkproPwNbgP9V4L1PrV3fMtIXiojjgdXAFOAF4L+klO4a6etKkjJHkCV1kv8MfC2l9HKB9z6VHGYfbMBrPQCcBJwGfBn4VkS8pQGvK0nCgCypTUXER2uHLPxhP88fHRHbI2JF7f4k4Czg0ibX+TcRkYBjyIdCvFKrO0XEbwz1cwCklHaklNallG5JKf0pcDvwB035QHvqGnLdktQuPMRCUrtaWbs+vZ/n/x4YD1xUu38WeXLe+tEurJdbgW+RzzaxCri67rlrgUNrtwf7OfoyDpg8oiqHbqj7X5LahgFZUru6FXiZfJjBXiLil4EK+Lu6U5/NB/7fgV40Ij4CTB9CHbenlH7U35MppcsjYjo5IH8rpXRZr/d7kiF8joj4G+BfgfXANOC9wDlAs8+FPNT9L0ltw9O8SWpbEXEdsBA4NKX0RO2x1wD3A5OAo1JKW4f4mo8Chw3hR76VUrpwgNf8CvDfgbemlNb28fygP0dEfBM4F3gDsBW4E/hcSumnQ6i5IUZj/0tSK/AYZEntrOdr/jPqHvskMBv4k+GEs5TSvJRSDOFy4SBe9hRgJ9DfmSYG/TlSShemlA5LKU1OKR2cUuoebDiOiEfrjn8ezGWgs300fP9LUivwEAtJ7awnoJ0G/DAijiFPVltNPu63uIiYABwP3JtS2t7PZs36HA8B24aw/c8HeL7l978kDYcBWVI7WwUk9kwU+xJ5YtjvpmEeP9boY5CB48jnK751H9s0/HP0JaW0pFGvVdOUuiWp2QzIktpWSunZiLgPOLW2Qt0S4MsppdtG8LIfYYjHIAP7Csgn1a77rWmUPseoa9e6JWkgBmRJ7e4G8ijtV4FngE+M5MVSSvMaUFO9g2rXzw2wXUM/RxO1a92S1C8n6Ulqdz3HwU4F/jSl9GzJYvrQs7T0ZyLicxHxqdpp0Hpr9c/Rn3atW5L6ZUCW1O4eqV2vAb5espC+pJRWAB8GXqxd/wV7Druo19KfYx/atW5J6pfnQZbU1iLiSvIiGaenlNaUrme42vVztGvdkrQvjiBLalu1iWG/SJ4Y1rbhrF0/R7vWLUkDcQRZUluJiLnk5ZWPAH4T+BkwP6X0UtHChqhdP0e71i1JQ+EIslpaRHwzIj5du312RDxQqI5PDWJVMTXH24C/Bn4J+P+A/9Sm4axdP0e71i1Jg2ZA1ojVlq99OSJeiIinaqF2aqPfJ6V0fUrp6EHUc2FE3NDo92+EiHhdRPyfiHgxIh6rfUXd37Yfi4i7I+L5iHgkIj7WzFpbVUrpstoSzwemlH4lpTTQam8tqV0/R7vWraGLiPdGxNpab38iIn4SEWcVrunRiOguWYM6gwFZjfKLKaWpwClAF3Bx7w1qS+62vFGu81JgB/B64NeAL0fEm/srhfwV9oHkUbuLIuKCUaxNkgCIiI8CXwT+ityv5gL/Ezh/iK/zqn5a6m9Bu/wNUmswIKuhUkobgZ8AbwGIiBQRvxsRPyMfq0hEvDMibo+If4+IVRFxQs/PR8TJEXFrbdT0++QlenueOyciNtTdnxMRP4yITRGxOSK+FBHHAl8BzqiNevx7bdsDIuKfats+FhEXR8S42nMXRsTKiPhCRGwGPtXPx5tUe43nI+KeiOgayr6JiNcA7wb+PKX0QkrpBuBK4Df62ZefTSndmlLalVJ6gPx19oKhvKckDVVEHABcQl4y/IcppRdTSjtTSj9OKX0sIiZHxBcj4ue1yxcjYnLtZ8+JiA0R8ScR8STwj7VD1K6IiP8VEc8BF9Z68tdrI9MbI+LTETG+robfiYj7av323og4JSK+TQ7qP6719z+ubXterSf/e0RcW/s70PM6j9ZquRN40ZCswTIgq6EiYg7wDvZeVvddwGnAcRFxMvAN4L+TVxj7KnBlreFOIi/Z+23gdcAPyIGyr/cZD/wL8BgwD5gFfC+ldB/wQWB1SmlqSml67Uf+HjgAeCOwiDwy+1t1L3ka8DB5pOQz/Xy884DvAdPJwfZLdfX8S60593X5l9pmRwG7UkoP1r3mHUB/I8j1nzeAs4F7BtpWkkboDPLgxP/p5/lPAKeTz+d9IjCfvb81fAO5hx8GfKD22PnAFeT++b+BbwK7gDcBJwP/CfhvAJEX0vkUuU+/ltx7N6eUfgN4nNo3limlz0bEUcB3yUvEzwSuIgfoSXX1vId8KsLpKaVdQ94b6kgGZDXKj2qjtTcA15G/luvx1ymlLSmll8nN8qsppZtSSrtTSt8CtpOb7enAROCLtdGKK8iLD/RlPnAo8LHa6Ma22ojsq9TC9AXkVb6eTyk9CnyevUduf55S+vvaaO3L/bznDSmlq1JKu8kh/sSeJ1JK70wpTe/n8s7aZlN59XLDW4Fp/bxfvU+R/73+4yC2laSROAh4Zh9h8teAS1JKT6eUNgF/yd799BXgL1JK2+v66eqU0o9SSq+QQ+87gI/U+vfTwBfIfRpyUP5sSmlNytallB7rp5ZfBf41pXR1Smkn8P8A+wFn1m3zdyml9fvo7dKr+FWDGuVdKaWl/Ty3vu72YcD7IuLDdY9NIofdBGxMe597sL+mOAd4bJCjATPIwbv+tR4jjzr3VWN/nqy7/RIwJSImDGFE4gXyH4Z6rwWe39cPRcRF5JGUs1NK2wf5XpI0XJuBGfvob4fy6n56aN39TSmlbb1+pvffgYnAE/nLMSAPAPRsMwd4aJC17lVLSumViFjP0Pu7tBdHkNUM9YF3PfCZXiOs+6eUvgs8AcyKuo5JPt6sL+uBuf0cT9b75N7PADvJTbn+dTfu42eGpDa7+4V+Lj+pbfYgMCEijqz70RPZx2ETEfHbwMeBJSmlDf1tJ0kNtJr8zd67+nn+57y6n9afzaSvftr778B2YEbd34HXppTeXPf8Ef28d+/X3quW2t+POTSwv6szGZDVbF8DPhgRp0X2moj4hYiYRm7Ku4Dfi4iJEfFfyYdS9OVmcqD+m9prTImInglsTwGze45Bqx0ScTnwmYiYFhGHAR8FGnZe45TS22vHxPV1eXttmxeBHwKX1GpeQD4u79t9vWZE/Br5UJUqpfRwo2qVpH1JKW0FPglcGhHvioj9az357RHxWfIxvxdHxMyImFHbdtD9NKX0BPBvwOcj4rURMS4ijoiIRbVN/gH4o4g4tfZ34k21vg25v7+x7uUuB34hIpZExETgD8nhe9UIdoFkQFZzpZTWAr9DnuD2LLAOuLD23A7gv9bubyEfW/bDfl5nN3mJ2zeRJ21sqG0PsJw8KvtkRDxTe+zDwIvkiXg3AN8hTxZstg+Rj497mvxH5n+klO6B/1gI5YW6bT9NPhZwTd1o9FeaXrGkjpNS+jx5IOFiYBN5VPci8kTqTwNrgTuBu4Bba48NxW+SD6+7l/y34ArgkNp7/4A8Wfo75EPQfkSe9Ad5kZqLaxOg/6h2hp9fJ0/Efob8d+EXa39PpGFzqWlJkiSpjiPIkiRJUh0DsiRJklTHgCxJkiTVMSBLkiRJdYotFDJjxow0b968Um8vSUXdcsstz6SUZpauw14sqZP114uLBeR58+axdu3aUm8vSUVFRH+rRDaVvVhSJ+uvF3uIhSRJklTHgCxJkiTVMSBLkiRJdQzIkiRJUh0DsiRJklRnwIAcEd+IiKcj4u5+no+I+LuIWBcRd0bEKY0vU5I6m71YkppnMCPI3wTeto/n3w4cWbt8APjyyMvah2eegeeeG9W3kKQW9E1apBenBI88MlqvLknlDRiQU0orgC372OR84J9SdiMwPSIOaVSBe3n4YZg5E77//VF5eUlqVa3Uiz/xCTjuONi2bTReXZLKa8QxyLOA9XX3N9Qee5WI+EBErI2ItZs2bRr6Ox1+OBx6KFx99bAKlaQxrGm9eMGCHI5vuGF4hUpSq2vqJL2U0mUppa6UUtfMmcNYYTUCqgqWLYPduxtfoCR1gJH24kWLYOJEWLp0FIqTpBbQiIC8EZhTd3927bHRUVWwZQvcfvuovYUktaGm9eKpU+GMM/wyT9LY1YiAfCXwm7UZ1KcDW1NKTzTgdfvW3Z2v7cySVK+pvbiq4Lbb8rxpSRprBnOat+8Cq4GjI2JDRLw/Ij4YER+sbXIV8DCwDvga8KFRqxbg9a+HE04wIEvqKK3Wi6sqn81i2bLRfBdJKmPCQBuklN4zwPMJ+N2GVTQY3d3wpS/BSy/B/vs39a0lqYRW68WnngoHHJDHKn71V5v1rpLUHO25kl5VwY4dcP31pSuRpI40YQIsXpwDckqlq5GkxmrPgLxwIUya5GEWklRQVcHjj8O6daUrkaTGas+AvP/++UScnmNIkoqpqnztWIWksaY9AzLkznzHHfDUU6UrkaSOdMQRMG+eAVnS2NPeARmcQi1JhfSs3bR8OezaVboaSWqc9g3IJ58MBx7o0IUkFdTdDc89B2vWlK5EkhqnfQPy+PGwZIlTqCWpoCVL8kiyYxWSxpL2DciQv9vbuBHuv790JZLUkQ46CE45xTnTksaW9g/IYGeWpIKqClavhuefL12JJDVGewfkww/P06j9bk+SiqmqPEnvuutKVyJJjdHeARlyZ772Wti5s3QlktSRzjwT9tvPsQpJY0f7B+Tu7vy93k03la5EkjrSlClw9tkGZEljR/sH5MWLYdw4O7MkFVRVcN99sGFD6UokaeTaPyAfeCB0dRmQJakg126SNJa0f0CG3Jlvvhm2bi1diSR1pOOPh4MPdqxC0tgwdgLy7t15sp4kqenGjctTQpYude0mSe1vbATkM86A/fd36EKSCqoqeOopuOuu0pVI0siMjYA8aRIsWmRAlqSCurvzta1YUrsbGwEZ8tDFgw/C44+XrkSSOtLs2XDMMQZkSe1vbAVksDNLUkFVBStWwPbtpSuRpOEbOwH5zW+GQw4xIEtSQVUFL78Mq1aVrkSShm/sBOSIfADcsmXwyiulq5GkjnTOOTBhgmMVktrb2AnIkIcunnkG7rijdCWS1JGmTYPTTzcgS2pvYysgO4VakoqrKrjlFti8uXQlkjQ8YysgH3JIPhbZgCxJxXR358VCli8vXYkkDc/YCsiQhy6uvz7PEpEkNd38+fDa1+ZV9SSpHY3NgLx9O9xwQ+lKJKkjTZgA557rl3mS2tfYC8iLFsHEiQ5dSFJBVQWPPAIPPVS6EkkaurEXkF/zGjjzTIcuJKkg126S1M7GXkCG3Jlvuw02bSpdiSR1pCOPhLlzDciS2tOgAnJEvC0iHoiIdRHx8T6enxsR10TEbRFxZ0S8o/GlDkHP6d6WLStahiQ1Ujv14p61m5Yvh927S1UhScMzYECOiPHApcDbgeOA90TEcb02uxi4PKV0MnAB8D8bXeiQdHXB9OkOXUgaM9qxF1cV/Pu/53MiS1I7GcwI8nxgXUrp4ZTSDuB7wPm9tknAa2u3DwB+3rgSh2H8eFi8OAfklIqWIkkN0na9eMmSfO1YhaR2M5iAPAtYX3d/Q+2xep8Cfj0iNgBXAR/u64Ui4gMRsTYi1m4a7eODqwrWr4ef/Wx030eSmqPtevHMmXDyyQZkSe2nUZP03gN8M6U0G3gH8O2IeNVrp5QuSyl1pZS6Zs6c2aC37odTqCV1npbrxVUFq1bBCy+M6ttIUkMNJiBvBObU3Z9de6ze+4HLAVJKq4EpwIxGFDhsRxwBhx9uQJY0VrRlL64q2LkTVqwoWYUkDc1gAvIa4MiIODwiJpEnflzZa5vHgSUAEXEsuSmXP8dadzdccw3s2lW6EkkaqbbsxQsWwOTJjlVIai8DBuSU0i7gIuCnwH3kGdL3RMQlEXFebbM/BH4nIu4AvgtcmFILzI6rKnjuObj55tKVSNKItGsv3m8/OPtsFzeV1F4mDGajlNJV5Akf9Y99su72vcCCxpbWAIsX55NxXn11Xl1PktpYu/biqoI/+RN44gk45JDS1UjSwMbmSno9DjoITj3V7/YkqaCeOdOOIktqF2M7IEPuzDfemA+1kCQ13Ykn5lO+OVYhqV10RkDevRuuu650JZLUkcaNy4uGLF3q2k2S2sPYD8hnnplniTh0IUnFdHfnY5Dvuad0JZI0sLEfkCdPhoULDciSVJDHIUtqJ2M/IEPuzPffDxs2lK5EkjrS3Llw1FGOVUhqD50TkMHOLEkFVVWeDrJjR+lKJGnfOiMgH388vP71frcnSQVVFbz4IqxeXboSSdq3zgjIEXmGyNKl8MorpauRpI50zjkwfrxf5klqfZ0RkCEH5KefhrvuKl2JJHWkAw6A+fMNyJJaX+cEZI9DlqTiqgrWroVnny1diST1r3MC8qxZcOyxBmRJKqiq8pFu11xTuhJJ6l/nBGTInXnFCti2rXQlktSRTjsNpk1zrEJSa+u8gLxtG6xaVboSSepIEyfmyXoGZEmtrLMC8qJFMGGCnVmSCqoqeOgheOSR0pVIUt86KyBPmwann25AlqSCurvztaeml9SqOisgQx66uPVW2Ly5dCWS1JGOOSbPm3asQlKr6syAnBIsW1a6EknqSBG5FS9bBrt3l65Gkl6t8wLyW9+az1bv0IUkFVNVsGUL3HZb6Uok6dU6LyBPmADnnpsDckqlq5GkjtRzHLJjFZJaUecFZMhDF489lqdRS5Ka7uCD4cQTDciSWlPnBmSwM0tSQVUFK1fCSy+VrkSS9taZAflNb4K5cw3IklRQdzfs2AHXX1+6EknaW2cG5J4p1MuXw65dpauRpI509tkwaZJjFZJaT2cGZMgBeetWWLu2dCWS1JH23x/OOsuALKn1dG5AXrIkjyTbmSWpmKqCO++Ep54qXYkk7dG5AXnGDDj5ZNc6laSCeuZM24oltZLODciQO/Pq1fDCC6UrkaSOdPLJ8LrXGZAltZbODsjd3bBzJ1x3XelKJKkjjRuXj3hz7SZJraSzA/JZZ8GUKR6HLEkFVRVs3Aj331+6EknKBhWQI+JtEfFARKyLiI/3s82vRMS9EXFPRHynsWWOkilT8nmGDMiSWtyY7cO4dpOk1jNgQI6I8cClwNuB44D3RMRxvbY5EvhTYEFK6c3AR0ah1tFRVXDvvXn4QpJa0Fjvw/Pm5fWbDMiSWsVgRpDnA+tSSg+nlHYA3wPO77XN7wCXppSeBUgpPd3YMkdRz9DFsmVl65Ck/o3tPkxuxddem6eFSFJpgwnIs4D1dfc31B6rdxRwVESsjIgbI+JtjSpw1J1wAsyc6dCFpFY2tvswOSC/8ALcdFPpSiSpcZP0JgBHAucA7wG+FhHTe28UER+IiLURsXbTpk0NeusR6plCvXSpU6gltbNB9WFozV587rm5HTtWIakVDCYgbwTm1N2fXXus3gbgypTSzpTSI8CD5Ea9l5TSZSmlrpRS18yZM4dbc+NVFTz5JNx9d+lKJKkvDevD0Jq9ePp0eOtbDciSWsNgAvIa4MiIODwiJgEXAFf22uZH5FELImIG+au+hxtY5+hyCrWk1jb2+zC5Fd98M2zdWroSSZ1uwICcUtoFXAT8FLgPuDyldE9EXBIR59U2+ymwOSLuBa4BPpZS2jxaRTfcnDlw9NEGZEktqSP6MDkg794N11xTuhJJnS5SoeNuu7q60tq1a4u8d58+/GH4xjdgyxaYPLl0NZLGuIi4JaXUVbqOVurFO3bkZaff9z649NLS1UjqBP314s5eSa9eVcFLL8Hq1aUrkaSONGkSnHNOnjMtSSUZkHssWgTjx3uYhSQV1N0NDz4Ijz9euhJJncyA3OOAA+C00wzIklSQc6YltQIDcr2qgrVr83HIkqSmO+44OPRQA7KksgzI9aoqLxayfHnpSiSpI0XkwyyWLYNXXildjaROZUCuN38+TJvm0IUkFVRV8MwzcPvtpSuR1KkMyPUmTszrnTqFWpKK6e7O17ZiSaUYkHvr7oaHH84XSVLTveEN8Ja3+GWepHIMyL05hVqSiqsquP56ePnl0pVI6kQG5N6OPhpmzzYgS1JBVQXbt8MNN5SuRFInMiD3FpE78/LlsHt36WokqSMtXJhX1nOsQlIJBuS+VBU8+yzcckvpSiSpI73mNXDmmQZkSWUYkPuyZEm+dgq1JBVTVflUb5s2la5EUqcxIPfl4IPhpJMcupCkgnrmTC9bVrYOSZ3HgNyf7m5YuRJefLF0JZLUkU45BQ480LEKSc1nQO5PVcHOnbBiRelKJKkjjR8PixfngJxS6WokdRIDcn/OPhsmT3boQpIKqipYvx4efLB0JZI6iQG5P/vtB2edZUCWpIJcu0lSCQbkfakquPtuePLJ0pVIUkd64xvzxZMKSWomA/K+9Axd2JklqZjubrjmGti1q3QlkjqFAXlfTjoJDjrI7/YkqaCqgueeg5tvLl2JpE5hQN6XcePyoiFOoZakYhYvhgjHKiQ1jwF5IFUFTzwB995buhJJ6kivex10dRmQJTWPAXkgTqGWpOKqCm68MR9qIUmjzYA8kMMOgyOPdKKeJBVUVbB7N1x7belKJHUCA/JgVFXuyjt2lK5EkjrSGWfA/vs7ViGpOQzIg9HdDS++mL/fkyQ13eTJsHChR7tJag4D8mCce24+o4WdWZKKqSq4/37YsKF0JZLGOgPyYEyfDvPnG5AlqSDnTEtqFgPyYFUVrFkDzz5buhJJ6khveQu84Q0GZEmjz4A8WFUFr7ziFGpJKiQiTwlZujS3Y0kaLYMKyBHxtoh4ICLWRcTH97HduyMiRURX40psEaefDlOnOnQhqRh7cQ7ImzbBXXeVrkTSWDZgQI6I8cClwNuB44D3RMRxfWw3Dfh94KZGF9kSJk6ERYsMyJKKsBdn3d352lYsaTQNZgR5PrAupfRwSmkH8D3g/D62+7+A/xvY1sD6WktVwbp18OijpSuR1HnsxcCsWXDccQZkSaNrMAF5FrC+7v6G2mP/ISJOAeaklP51Xy8UER+IiLURsXbTpk1DLrY4p1BLKsdeXFNVsGIFbBuT/wsgqRWMeJJeRIwD/hb4w4G2TSldllLqSil1zZw5c6Rv3XzHHguHHmpAltRyOqkXV1UOxytXlq5E0lg1mIC8EZhTd3927bEe04C3ANdGxKPA6cCVY3FyCBG5My9b5hRqSc1mL65ZtAgmTHCsQtLoGUxAXgMcGRGHR8Qk4ALgyp4nU0pbU0ozUkrzUkrzgBuB81JKa0el4tKqCrZsgdtuK12JpM5iL66ZOhXOOCOf7k2SRsOAATmltAu4CPgpcB9weUrpnoi4JCLOG+0CW86SJfnaoQtJTWQv3ltVwa23wubNpSuRNBZFSqnIG3d1daW1a9t0YOOEE2DmzHyohSQNQ0TcklIqfvhDu/biG2/Mo8jf/z78yq+UrkZSu+qvF7uS3nBUFdxwA7z0UulKJKkjdXXBAQf4ZZ6k0WFAHo6qgh074PrrS1ciSR1pwgRYvDgH5EJfhEoawwzIw7FwIUya5AwRSSqoquCxx/L6TZLUSAbk4dh/f1iwwO/2JKkgl52WNFoMyMNVVXDHHfDUU6UrkaSO9KY3wWGH+WWepMYzIA9Xz9CFZ7KQpCJ61m5avhx27SpdjaSxxIA8XKecAgce6Hd7klRQVcHWrdCGZ6qT1MIMyMM1fnxeNMQp1JJUzJIleSTZsQpJjWRAHomqgo0b4f77S1ciSR3poIPyF3oGZEmNZEAeiarK184QkaRiurth9Wp4/vnSlUgaKwzII3H44XDEEQ5dSFJBVZUn6a1YUboSSWOFAXmkurvh2mth587SlUhSR1qwAKZMcaxCUuMYkEeqqvL3ejfdVLoSSepIU6bkBU4NyJIaxYA8UosXw7hxdmZJKqiq4N5787xpSRopA/JIHXggdHUZkCWpoJ61m5wzLakRDMiNUFVw8835bPWSpKY74QSYOdOxCkmNYUBuhKqC3bvzZD1JUtONG5dHkZcude0mSSNnQG6E00+H/fd36EKSCqoqeOopuPvu0pVIancG5EaYPBkWLTIgS1JBPWs32YoljZQBuVGqCh58EB5/vHQlktSRZs+GY44xIEsaOQNyozh0IUnFdXfDddfB9u2lK5HUzgzIjfLmN8Mhh3iOIUkqqKrg5Zdh1arSlUhqZwbkRonYM4X6lVdKVyNJHemcc2D8eL/MkzQyBuRG6u6GZ56BO+4oXYkkdaTXvjafWMgv8ySNhAG5kXqWcnLoQpKKqSpYuxa2bCldiaR2ZUBupEMPzcciG5AlqZiqyouFLF9euhJJ7cqA3GhVBddfn2eJSJKa7q1vhWnTHKuQNHwG5Earqnx+oZUrS1ciSR1p4kQ491wDsqThMyA32sKFuTvbmSWpmKqCRx6Bhx4qXYmkdmRAbrSpU+GMMwzIklRQz9pNns1C0nAMKiBHxNsi4oGIWBcRH+/j+Y9GxL0RcWdELIuIwxpfahupKrjtNti0qXQlksYI+/DQHHUUzJnjWIWk4RkwIEfEeOBS4O3AccB7IuK4XpvdBnSllE4ArgA+2+hC20rP0MWyZWXrkDQm2IeHLiK34uXLYffu0tVIajeDGUGeD6xLKT2cUtoBfA84v36DlNI1KaWXandvBGY3tsw20/4PpnkAAA2ESURBVNUF06c7dCGpUezDw9DdDc8+C7fcUroSSe1mMAF5FrC+7v6G2mP9eT/wk76eiIgPRMTaiFi7aSwffjB+PCxenANySqWrkdT+GtaHoXN68ZIl+dqxCklD1dBJehHx60AX8Lm+nk8pXZZS6kopdc2cObORb916qgrWr4ef/ax0JZI6yEB9GDqnFx98MJx0kgFZ0tANJiBvBObU3Z9de2wvEdENfAI4L6W0vTHltTGXnZbUOPbhYaoqWLUKXnihdCWS2slgAvIa4MiIODwiJgEXAFfWbxARJwNfJTflpxtfZhs64giYN8+ALKkR7MPDVFWwc2de4FSSBmvAgJxS2gVcBPwUuA+4PKV0T0RcEhHn1Tb7HDAV+EFE3B4RV/bzcp2jZwr1NdfArl2lq5HUxuzDw3fWWTB5smMVkoZmwmA2SildBVzV67FP1t3ubnBdY0NVwde+BjffDGeeWboaSW3MPjw8++2XQ7IBWdJQuJLeaFq8OI8ku5STJBVTVXD33fDEE6UrkdQuDMij6aCD4NRTHbqQpIJcdlrSUBmQR1tVwY03wvPPl65EkjrSSSfBjBmOVUgaPAPyaOvuzpP0rr22dCWS1JHGjcuLhixd6tpNkgbHgDzaFizIs0QcupCkYqoqH4N8772lK5HUDgzIo23yZFi40IAsSQW5dpOkoTAgN0NVwf33w4YNpSuRpI502GFw5JEGZEmDY0BuBqdQS1JxVQXXXQc7dpSuRFKrMyA3w/HHw8EHO3QhSQVVFbz4IqxeXboSSa3OgNwMEfkAuKVL4ZVXSlcjSR3p3HPzGS0cq5A0EANys1QVPP003HVX6UokqSMdcADMn+/RbpIGZkBulp7jkB26kKRiqgrWrIFnny1diaRWZkBullmz4NhjHbqQpIKqKh/pds01pSuR1MoMyM1UVbBiBWzbVroSSepIp58OU6f6ZZ6kfTMgN1NVwcsvw6pVpSuRpI40cSKcc44BWdK+GZCbadEimDDBzixJBVUVPPQQPPJI6UoktSoDcjNNm5a/3zMgS1IxPctOOyVEUn8MyM1WVXDrrbB5c+lKJKkjHXssHHqoYxWS+mdAbraqgpRg+fLSlUhSR4rIrXjZMti9u3Q1klqRAbnZ3vrWfLZ6hy4kqZiqgi1b4LbbSlciqRUZkJttwoS83unVV+eRZElS0/Uch+xYhaS+GJBL6O6GRx/N06glSU33+tfDCScYkCX1zYBcgstOS1Jx3d2wciW89FLpSiS1GgNyCUceCXPnNj8gb9sGP/95c99TklpUVcGOHXD99aUrkdRqDMgl9EyhXr589KdQ796dT/b5W7+Vv1OcNQtOOQX++q9h3brRfW9JamELF8KkSX6ZJ+nVDMilVBVs3Qpr1zb+tVOCNWvgIx+B2bPze/3wh/Dud+dgPGkS/Nmf5ZFsw7KkDrX//rBggQFZ0qsZkEtZsiRfN7IzP/ggfOpTcPTRMH8+fPnLcMYZcMUV8OST8I1vwMc/DjfeCI89Bp//PEyevCcsn3wy/NVfwc9+1riaJKmFVRXceSc89VTpSiS1EgNyKTNm5EA60oD8xBPwhS/k8ysffTRccgnMmQP/8A+54/eMHO+3394/N3cufPSjsHp1Dst/+7cwZQp84hNw1FGGZUkdoWfOtMtOS6pnQC6pqnJAfeGFof3c1q3wj/+Yp2DPnp2Dbkp5RHjDhrw81PvfD9OnD+715s6FP/iDXMvjj+ewvN9+e8LySSfBZz6TR6glaQw5+WQ48EAPs5C0NwNySVUFO3fCddcNvO22bXtGg1//evjt387nUr74Yrjvvnws80c/CoceOrKa5szJYXnVqhyWv/CFfKDexRfnEWrDsqQxZPz4fMTb0qWu3SRpDwNySWedlQ9r6G/oYvfuPaPBb3hDDscrV8IHPwg33ZQPf/jLv4Rjjhmd+ubMyRP96sPya16zJyyfeCJ8+tPwwAOj8/6S1ARVBRs3wv33l65EUqsYVECOiLdFxAMRsS4iPt7H85Mj4vu152+KiHmNLnRMmjIFzj5774PfUoJbbsmjwXPm5MMofvADeNe74N/+LR9C8cUv5kl4Ec2rtScsr1wJ69fnGqZOhT//8xzQDcvSqLMXjw7XbpLUW6QBvlOKiPHAg0AFbADWAO9JKd1bt82HgBNSSh+MiAuA/5JS+tV9vW5XV1daOxqnOGs3n/sc/PEfw4oVcM018J3v5JA5aRK84x3w3vfCO9/56kl2rWLDBvjnf84hfuXK/Njxx8Mv/zIcfvje2/YO9Pu6P5Rte+73dxno+aFu23ubRt/v/Vj9523WY7338UDbDPd6qO83ku1bTETcklLqGsL29uJR9KY3wbHHwo9//OrnUoLt2/Nl27bGXL/yCowbt+efYM/t/q5H8lz9pb/Hh7tdX9tC49rtQNv2fnyg+8P5mZ779Y+P9mPNuB7JYyN5jVbTXy+eMIifnQ+sSyk9XHuh7wHnA/fWbXM+8Kna7SuAL0VEpIHSt/YMXSxcmH97Fi2CP/qjfDjFgQeWrW0wZs+G3//9fNm4MYflyy+HT36ydGXSvvXu1v118oFu3357Tlajz148iqoKvv71PC+5d5jdvr0x7zFpUj6z5uTJOVSmlINySnvf7n3d32PSWDDc1lt/e8GCvPZaIw0mIM8C1tfd3wCc1t82KaVdEbEVOAh4pn6jiPgA8AGAuXPnDrPkMeaEE+BjH4ODD4YLLsiBs13NmgW/93v5smlTPttGj97dfF/3h7Jtz/3+LgM9P9Rte2/T6Pu9H6v/vM16rPc+Hmib4V4P9f1Gsv2+bg/nZ+pvH3QQTWIvHkUf+hBs3pwn7U2ZkkPsvq4Hs039tpMm5VDcSEMJ2PWX/h4f7nb12/bU1Yh2O9C2vR8f6P5IWnT946P9WDOuR/LYSF5jtG4fdhgNN5iA3DAppcuAyyB/rdfM925Z48bBZz9buorGmzkzXyS1HHvxqx1/fP7yq51E5EAvqfEG8/+zG4E5dfdn1x7rc5uImAAcAGxuRIGSJMBeLElNM5iAvAY4MiIOj4hJwAXAlb22uRJ4X+32LwHLPeZNkhrKXixJTTLgIRa149guAn4KjAe+kVK6JyIuAdamlK4Evg58OyLWAVvIjVuS1CD2YklqnkEdg5xSugq4qtdjn6y7vQ345caWJkmqZy+WpOZwJT1JkiSpjgFZkiRJqmNAliRJkuoYkCVJkqQ6UeoMQBGxCXhsmD8+g14rQ3Uo90PmftjDfZG1w344LKVUfDWdEfTidtjHzeK+yNwPmfsha5f90GcvLhaQRyIi1qaUukrXUZr7IXM/7OG+yNwPo899vIf7InM/ZO6HrN33g4dYSJIkSXUMyJIkSVKddg3Il5UuoEW4HzL3wx7ui8z9MPrcx3u4LzL3Q+Z+yNp6P7TlMciSJEnSaGnXEWRJkiRpVBiQJUmSpDotG5Aj4m0R8UBErIuIj/fx/OSI+H7t+ZsiYl7zq2yOQeyLCyNiU0TcXrv8txJ1jqaI+EZEPB0Rd/fzfETE39X20Z0RcUqza2yWQeyLcyJia93vwyebXWMzRMSciLgmIu6NiHsi4vf72KZjfi9Gi704sw9n9uLMPpyN6T6cUmq5CzAeeAh4IzAJuAM4rtc2HwK+Urt9AfD90nUX3BcXAl8qXeso74eFwCnA3f08/w7gJ0AApwM3la654L44B/iX0nU2YT8cApxSuz0NeLCPfxsd83sxSvvYXjz4/TDm+3Dtc9qLB7cf7MNt/jvRqiPI84F1KaWHU0o7gO8B5/fa5nzgW7XbVwBLIiKaWGOzDGZfjHkppRXAln1scj7wTym7EZgeEYc0p7rmGsS+6AgppSdSSrfWbj8P3AfM6rVZx/xejBJ7cWYfrrEXZ/bhbCz34VYNyLOA9XX3N/DqHf4f26SUdgFbgYOaUl1zDWZfALy79tXFFRExpzmltZTB7qdOcUZE3BERP4mIN5cuZrTVvtY/Gbip11P+XoyMvTizDw+e/+b2sA9nbfk70aoBWUPzY2BeSukE4Gr2jOaoM91KXlv+RODvgR8VrmdURcRU4J+Bj6SUnitdjzqWfVj17MNtrlUD8kag/v++Z9ce63ObiJgAHABsbkp1zTXgvkgpbU4pba/d/Qfg1CbV1koG8zvTEVJKz6WUXqjdvgqYGBEzCpc1KiJiIrkp/++U0g/72MTfi5GxF2f24cHz3xz24V7a8neiVQPyGuDIiDg8IiaRJ35c2WubK4H31W7/ErA81Y4GH2MG3Be9juU5j3wMUKe5EvjN2mzZ04GtKaUnShdVQkS8oecY0IiYT/53PtYCC7XP+HXgvpTS3/azmb8XI2MvzuzDg+e/OezDvbTl78SE0gX0JaW0KyIuAn5Knj38jZTSPRFxCbA2pXQl+T/ItyNiHflA+QvKVTx6Brkvfi8izgN2kffFhcUKHiUR8V3yrOAZEbEB+AtgIkBK6SvAVeSZsuuAl4DfKlPp6BvEvvgl4H9ExC7gZeCCMRhYABYAvwHcFRG31x77M2AudN7vxWiwF2f24T3sxZl9+D+M2T7sUtOSJElSnVY9xEKSJEkqwoAsSZIk1TEgS5IkSXUMyJIkSVIdA7IkSZJUx4AsSZIk1TEgS5IkSXX+f3ppEbsYXmd3AAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 720x288 with 2 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "collapsed": true,
        "id": "rm1aCQ_KthEt"
      },
      "source": [
        ""
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}